home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
c
/
jpl_c.zip
/
RANDOM.C
< prev
next >
Wrap
Text File
|
1986-05-18
|
5KB
|
138 lines
/* 1.1 10-16-84 */
/************************************************************************
* Robert C. Tausworthe *
* Jet Propulsion Laboratory *
* Pasadena, CA 91009 1984 *
************************************************************************
* Uniform Random Number Generator, programmed using the algorithm
* given in:
*
* Tausworthe, Robert C., STANDARDIZED DEVELOPMENT OF COMPUTER
* SOFTWARE, Part 2, Appendix L, Example 3, Prentice-Hall, Inc.,
* Englewood Cliffs, NJ, 1979, pp. 499-506.
*
* The generator uses the primitive polynomial over GF(2^521),
* f(x) = x^521 + x^32 + 1, found in Zierler, N., and Brillhart, J.,
* "On Primitive Polynomials (Mod 2), II," INFORMATION AND CONTROL,
* Volume 14, No. 6, pp. 566-569, June, 1969.
*
* Numbers output by the generator are unsigned, uniformly distributed
* over 0..65535. Numbers have no serial statistical correlation,
* adjacent numbers are uniformly distributed over the 32-cube, and
* runs-up/down of length up to 13 require more than 10^11 samples to
* show deviation from theoretical. The period of the generator is
* slightly larger than 1.3 * 10^154.
*
* For more information on Tausworthe random number generators, see;
*
* Tausworthe, Robert C., "Random Numbers Generated by Linear Recurrence
* Modulo 2," MATH COMP., Vol. XIX, No. 90, pp. 201-208, April 1965.
*
* Toothill, J. P. R., "An Asymptotically Random Tausworthe Sequence,"
* ACM J., Vol. 20, No. 3, pp. 469-481, July 1973.
*
* Toothill, J. P. R., et al., "The Runs-Up and -Down Performance of
* Tausworthe Pseudorandom Number Generators," ACM J., Vol. 18, No. 3,
* pp. 381-399, July 1971.
*
*----------------------------------------------------------------------*
*
* The shift register of 521 bits requires 65.125 bytes, or 65 bytes
* plus one bit. One register load then provides 32 16-bit numbers
* before it must be refilled. The 32-bit tap is at a 4-byte shift
* from the first byte, and at 521 - 32 = 489 = (61 bytes + 1 bit)
* inverse shift from the end.
*
*----------------------------------------------------------------------*/
#include "defs.h"
#include "stdtyp.h"
/*----------------------------------------------------------------------*/
#define MULTIPLIER 18829 /* This is 5^15 (mod 65536) */
#define PRESEED 0xe5a1 /* a number picked out of the air! */
/*\p--------------------------------------------------------------------*/
LOCAL utiny shiftreg[66]; /* 521-bit shift register */
LOCAL utiny *randptr = shiftreg - 1; /* initialized for not accessed */
LOCAL utiny *endreg = shiftreg + 62; /* always points at 32nd number */
/************************************************************************/
VOID
randize(seed) /* Randomize the generator using seed.
See algorithm for seeding values. */
/*----------------------------------------------------------------------*/
{
int i, *rseed;
if (seed IS 0)
seed = PRESEED;
else if (seed < 0)
{ rseed = (int *) &shiftreg[16];
while ((seed = *(++rseed)) IS 0)
; /* find something non-zero, whatever. */
}
for (i = 0; i < 66; i++)
/* assume no adverse reaction on integer overflow */
{ shiftreg[i++] = ((seed *= MULTIPLIER) &0xff);
shiftreg[i] = ((seed >> 8) & 0xff);
}
shiftreg[65] &= 0x80; /* mask out bits beyond 521st. */
refill(); /* refill and reset randptr */
}
/************************************************************************/
unsigned
urandom() /* Return the next random number. If shift register is
depleted, refill the entire register. */
/*----------------------------------------------------------------------*/
{
unsigned r;
if (randptr < shiftreg)
randize(0); /* in case randize wasn't called first */
else if (randptr > endreg)
refill(); /* and reset randptr & randx */
r = UTINY(*randptr++);
r += (UTINY(*randptr++) << 8);
return (r);
}
/*\p*********************************************************************/
double
random() /* Return a uniform random value in [0, 1]. */
/*----------------------------------------------------------------------*/
{
unsigned urandom();
return (urandom() / 65535.0);
}
/************************************************************************/
LOCAL VOID
refill() /* refill the shift register */
/*----------------------------------------------------------------------*/
{
utiny *p, *q, cy0, cy1;
FAST int i;
p = shiftreg; /* point at 1st byte. */
q = &p[4]; /* shift 4 * 8 = 32, */
for (i = 0; i < 62; i++) /* and mod-2 add the register */
*p++ ^= *q++; /* to the 32-bit shift. */
p--; /* p should now be at byte 61: */
q = shiftreg; /* first 489 bits are ok. */
cy0 = 0; /* set carry bit to zero */
for (i = 0; i < 4; i++) /* for the remaining 32 bits: */
{ cy1 = ((*q & 1) ? 0x80 : 0); /* save carry for nxt shift*/
*p++ ^= (((*q++ >> 1) & 0x7f) | cy0);
/* mod-2 add the 489-shifts */
cy0 = cy1; /* set carry for next shift */
}
*p ^= cy0; /* and mod-2 the final bit. */
randptr = shiftreg; /* point at the first number */
}